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Structural defects in materials such as vacancies, grain boundaries, and dislocations may trap 
hydrogen and a local accumulation of hydrogen at these defects can lead to the degradation of the 
materials properties. An important aspect in obtaining insight into hydrogen induced embrittlement 
on the atomistic level is to understand the diffusion of hydrogen in these materials. In our study we 
employ kinetic Monte Carlo (kMC) simulations to investigate hydrogen diffusion in bcc iron within 
different microstructures. All input data to the kMC model, such as available sites, solution energies, 
and diffusion barriers are obtained from first-principles calculations. We find that hydrogen mainly 
diffuses within the interface region with an overall diffusivity that is lower than in pure bcc-Fe bulk. 
The concentration dependence of the diffusion coefficient is strongly non-linear and the diffusion 
coefficient may even decrease with increasing hydrogen concentration. To describe the macroscopic 
diffusion coefficient we derive an analytic expression as a function of hydrogen concentration and 
temperature which is in excellent agreement with our numerical results for idealised microstructures. 
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I. INTRODUCTION 

Hydrogen is a very common impurity in iron-based ma- 
terials. It is incorporated into the material during pro- 
duction and service and exhibits a high mobility within 
the bulk phase. It has further been shown that mi- 
crostructural defects in the material such as vacancies, 
dislocations, and grain boundaries can trap hydrogen 
impurities^™^. The resulting local accumulation of hy- 
drogen at these defects can then lead to a degradation 
of the mechanical properties of the material, which is 
also referred to as hydrogen embrittlementP. To ex- 
plain the complex mechanisms underlying hydrogen em- 
brittlement various approaches have been developed in- 
cluding hydrogen-enhanced local plasticity (HELP jPUl^ 
hydrogen-enhanced decohesion (HEDE) 9 "^ and super- 
abundant vacancy formation^!. 

The most important aspects in trying to understand 
the behaviour of hydrogen within a material and its role 
in hydrogen embrittlement are the hydrogen solubility, 
the interaction of hydrogen with defects, and the hydro- 
gen mobility. The solubility of hydrogen within pure bcc- 
Fe is relatively low. Point and extended defects, however, 
can provide interstitial sites that are energetically much 
more favourable for hydrogen than the tetrahedral site in 
bcc-Fe. In a recent study we investigated the solution en- 
ergy of H in the presence of open and close-packed grain 
boundaries (GB) in bcc- and f ee- Fe^ employing density- 
functional theory (DFTjpEU calculations. We find that 
in general hydrogen prefers the interstitial sites within 
the grain boundary region (with the exception of close- 
packed grain boundaries in fcc-Fe). In particular we find 



a large energy gain of 0.4 — 0.5 eV for hydrogen intersti- 
tials at the E5[001](310) grain boundary in bcc-Fe (£5 
in short) as compared to the bulk region. This indicates 
that hydrogen segregates to the grain boundary and is 
trapped there. 

The mobility of hydrogen within pure bcc-Fe bulk 
is high. The calculated diffusion barrier for H atoms 
movin g be tween neighbouring tetrahedral sites is ~ 
0.1 e V 15 * 18 l This is consistent with earlier experimen- 
tal observations indicating that H diffuses rapidly within 
bcc metalsSsHUl Experimental studies investigating the 
diffusion of hydrogen in the presence of dislocations, 
grain boundaries, and phase boundaries show that the 
measured diffusivity depends on the hydrogen concentra- 
tion. In bulk Pd, Pd grain boundaries, and Pd/Al 2 03 
phase boundaries, the diffusion constant generally in- 
creases with increasing H concentration^. Within a the- 
oretical study employing kinetic Monte Carlo (kMC) sim- 
ulations it was illustrated that the attraction of H to 
screw dislocations in Fe materials can significantly affect 
H diffusion 2 . To properly evaluate the mobility of hydro- 
gen within a certain material it is thus important to also 
consider the effect of the microstructure. 

In this study, we employ kMC simulations to address 
H diffusion in the presence of grain boundaries in bcc- 
Fe. In a previous study we found that the H diffusion 
barriers within the grain boundary region of bcc-Fe are 
much higher (0.25 — 0.6 eV)pSl than in the bulk region. 
This suggests that H interstitials diffuse relatively slowly 
or are effectively immobile within the grain boundary in- 
terface, and that therefore the grain boundaries do not 
provide fast diffusion channels for hydrogen. Utilising 
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the information about available interstitial sites, diffu- 
sion barriers, and solution energies extracted from our 
DFT calculations, we set up a series of kMC models that 
represent idealised microstructures. The first model is an 
idealised cubic grain structure in bcc-Fe, with and with- 
out point defects included in the bulk region of the grain. 
Since the model grains are much smaller than grains in 
the actual material we consider in our second model a 
parallel arrangement of grain boundary interfaces. This 
layered structure naturally introduces an anisotropy in 
the diffusivity. The third model presented in this paper 
represents more detailed the structure of the £5 grain 
boundary in bcc-Fe. Within the above models, H dif- 
fusion constants are determined as a function of H con- 
centrations and temperatures. From the results obtained 
within these model systems general trends for H diffusiv- 
ity in various structural environments can be extracted. 
We then compare the numerical results to the derived 
analytic expression that describes the macroscopic diffu- 
sivity of hydrogen in different microstructures. 

The paper is organised as follows. The computational 
approach is detailed in Sec. [ll| In Sec. |III| and |IV| the 
results for H diffusion within the idealised cubic grain 
structure with and without additional point defects are 
presented. The results for the layered structure and for 
the more detailed model of the £5 grain boundary are 
discussed in Sec. [V] an d |VI| respectively. Our findings 
are summarised in Sec. IVII1 



II. COMPUTATIONAL APPROACH 



We employ kinetic Monte CarlcE^U simulations to in- 
vestigate H diffusion under various conditions and within 
a number of idealised microstructures. Within kMC sim- 
ulations the time evolution of the system is described by a 
stochastic trajectory. The system states along this trajec- 
tory are connected by processes associated with a certain 
probability. Here, we use a lattice approach, i.e. possi- 
ble atomic positions are mapped onto a lattice and the 
system can evolve by atoms (hydrogen) hopping between 
neighbouring lattice sites (interstitial sites in bcc-Fe bulk 
and grain boundaries) . Within harmonic transition state 
theory^ the microscopic rate constant for a hop, ki, as- 
sociated with process i can be written as: 



h = i>o,« CX P (-A£7i/fcBT) 



(1) 



where vqj is the attempt frequency, &b is the Boltz- 
mann constant, T is the temperature, and AEi is the 
energy barrier associated with process i. For all pro- 
cesses that can occur within our kMC models the corre- 
sponding energy barriers were calculated employing DFT 
calculations (details regarding the DFT calculations can 
be found in Ref. I15|) . For a given system configuration 
the sum over all rate constants of all possible processes is 
evaluated, fctot = ^i, an d a process p to move to the 



next system state is chosen according to 
p-i p 



(2) 



i=0 



i=0 



where p\ is a uniform random number between and 1. 
Since the kMC algorithm simulates a sequence of Poisson 
processes, the real time evolution for each kMC step can 
be evaluated as^ 



t-tt- Mp 2 )/k toU 



(3) 



where p 2 is a second random number between and 1. 

The diffusion constant tensor is calculated from the 
mean square displacement of the hydrogen atoms. To 
obtain better statistics on the diffu sion co nstants we fol- 
low an approach outlined previousl y^ 24 * 29 !, and divide the 
kMC trajectory into a number of segments. The diffu- 
sion constant is calculated as the time weighted average 
of the diffusion constants for each segment i 



D kk = D kk ,iAti/t 



(4) 



where for the tetragonal structures analysed in this pa- 
per the diagonal components of the diffusion tensor Dkk 
suffice; Ati — ti — is the time length of segment i, t 
is the total length of the kMC trajectory, and 



(5) 



is the diffusion constant for segment i. Here, rk(U) is the 
position of a H atom in fc-direction (with k — x,y,z) at 
time ti and (...) denotes the average over all particles. 
The overall diffusion constant D is is defined as the av- 
erage of the diagonal components of the diffusion tensor, 

D=\T,kDkk. 

Similarly, we define the probability p s of finding hy- 
drogen atoms in certain sites s (e.g. bulk, interface or 
point defect sites) as 



J2iPi,sAU 



with 



N, 



Pi 



i,H,s 



H 



(6) 



(7) 



where TVh is the total number of hydrogen atoms in the 
simulation and iVi h,» is the number of hydrogen atoms 
at site type s in the trajectory segment i. 

The analytic expression for the diffusion constant D 
within a symmetric lattice and in the limit of low hy- 
drogen concentrations is given by the classical Arrhenius 
expression 



D = D cxp (-AE/k B T) 



(8) 



with Do = Ta^vo- AE is the diffusion barrier, ao is the 
jump distance for a diffusion hop, and T is a geomet- 
ric prefactor that is related to the connectivity of each 
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TABLE I: Diffusion barriers, AEi 
within the idealised cubic grain. 



between various sites 



process i 



bulk site — > bulk site 
bulk site — > interface site 
interface site — > bulk site 
interface site — ► interface site 
bulk site — > point defect 
point defect — > bulk site 



AE t (eV) 



0.088 
0.088 
0.600 
0.250 
0.088 
0.750 



FIG. 1: (Colour online) Model of an idealised cubic grain 
in bcc-Fe. The spheres represent tetrahedral interstitial sites 
in bcc-Fe, each having four nearest neighbours at the same 
distance. Interface sites are shown in blue, bulk sites are 
shown in green. 



site to its neighbouring sites. For isotropic diffusion the 
geometric prefactor follows^ 



r = 



2d 



(9) 



where n is the number of equivalent nearest neighbour 
sites and d is the dimensionality of the system. 

To obtain an analytic expression for the probability 
to find hydrogen in various types of interstitial sites we 
assume a grand canonical ensemble and exp ress the total 
number of H atoms within the system aa^ 



exp((£ s - M )/fc B T) 



(10) 



where the sum runs over the various interstitial site types 
s, n s is the number of sites of type s, E s is the solution 
energy of hydrogen at site type s, and \x is the chemical 
potential. The probability to find hydrogen at site type 
s is therefore 



Ps 



iV H 



1 



1 + cxp ((£ s - n)/kaT) 



(11) 



The analytic results for the probability of finding hydro- 
gen in certain sites can directly be compared with our 
numerical results employing Eq. Q. 



III. IDEALISED CUBIC GRAIN 

To study hydrogen diffusion within an idealised cu- 
bic grain structure of bcc-Fe we use the model shown in 
Fig. [T] In this simple model we only consider two dif- 
ferent site types, bulk and interfaces sites, shown as blue 
and green spheres in Fig. [T] respectively. All sites in the 
kMC model are arranged as tetrahedral sites in bcc-Fe 
(preferred interstitial site for hydrogen), each having four 
nearest neighbours at the same distance. Thus the two 
site types have the same geometry, but different solu- 
tion energies and respective diffusion barriers. The grain 
boundaries between bulk regions are represented by three 



layers of interface sites. Within the kMC simulations we 
use a simulation cell with a side length of x — 45.312 A 
containing 5952 interface and 43200 bulk sites. The simu- 
lation cell is periodically repeated in all three dimensions. 
The total number of H atoms within the system is 



H.intf 



H.bulk 



(12) 



where A^mtf and A?H,buik are the number of H atoms at 
the interface and bulk sites, respectively. Correspond- 
ingly, we express the number of H atoms within the two 
site types as 



"intf 



intf 



and 



1 + exp {(E int{ - n)/k B T) 



"bulk 



(13) 



H.bulk — 



1 + exp ((-Ebuik - fj)/k B T) 
The probability of finding H at bulk sites is therefore 



(14) 



N- 



Pbulk — 



H,bulk 



^H.bulk + iVH.intf 



(15) 



To determine the solution energies and respective diffu- 
sion barriers between the two site types we use the results 
of our recent DFT study^. The bulk sites in our model 
correspond to tetrahedral interstitial sites in bcc-Fe, and 
the interface sites correspond to the most stable intersti- 
tial site within the grain boundary region of the £5 GB. 
According to our DFT calculations the interstitial site 
at the £5 GB provides a trap for H interstitials with a 
binding energy of AE hind = ^buik - E inti = 0.512 eV. 
There are four different processes for hydrogen diffusion 
between bulk and interface sites: bulk — > bulk, bulk — > in- 
terface, interface — > bulk, and interface — > interface. The 
energy barriers, AEi, to determine the rate constants for 
these processes according to Eq. ([I]) are likewise taken 
from our recent DFT studjff^ and are listed in Table [i] 
The attempt frequency v$ is set to 10 13 s _1 for all pro- 
cesses. 

In Fig. [2] results for hydrogen diffusion constants for 
the cubic grain structure and for a perfect bulk structure 
are compared within the dilute limit. The perfect bulk 
structure contains only bulk sites with a diffusion bar- 
rier of 0.088 eV between neighbouring sites and serves 
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FIG. 2: (Color online) Logarithm of the diffusion constants 
for hydrogen as a function of the inverse temperature within 
the dilute limit. The black triangles are the results for hy- 
drogen diffusion between tetrahedral sites in perfect bcc-Fe. 
The blue, green, and red symbols represent the diagonal com- 
ponents of the diffusion tensor for the idealised cubic grain 
structure confirming an isotropic diffusion of hydrogen. The 
dashed lines correspond to linear fits of the kMC data. 



as a reference within our simulations. In Fig. [2] the log- 
arithm of the diffusion constants is shown as a function 
of the inverse temperature. If we assume an Arrhenius- 
like behaviour, cf. Eq. ([8]), then the slope corresponds to 
—AE/kB and the y-intercept to ln(D ). For our refer- 
ence, perfect bulk system a barrier of AE = 0.088 eV is 
extracted from the slope. We obtain ln(Dn) = —16.52, 
as expected from the theoretical value of ln(raQi'o) with 
r = 4/6, a = 1.001 A, and u = 10 13 . Within the grain 
structure the three diagonal components of the diffusion 
tensor, D xx , D yy , and D ZZl are equivalent confirming the 
isotropic diffusion of hydrogen. The slope of the linear 
fit to the simulation data yields a value for the effective 
diffusion barrier of AE — 0.256 eV, which corresponds to 
the barrier for an interface — > interface hop, suggesting 
that diffusion mainly takes place within the interface re- 
gion. The fitted value of ln(Dn) = —16.95 is smaller than 
in perfect bulk. Assuming diffusion only within the inter- 
face region there are only two adjacent interface sites, i.e. 
r = 2/6, yielding a theoretical value of ln(Z?n) = —17.21. 
It can also be seen that within the dilute limit diffusion 
is slower in the grain structure than in the perfect bulk 
structure. This can be explained by the fact that within 
the grain structure hydrogen is confined to the interface 
region which exhibits a higher barrier for diffusion. 

In Fig. [3] hydrogen diffusion constants as a function of 
H concentration are shown for temperatures of T = 300, 
600, and 1000 K, respectively. The kMC results are illus- 
trated by blue triangles. For all three temperatures the 
diffusivity is constant for low H concentrations, it then 
decreases as the number of H atoms approaches the num- 
ber of interface sites, before it again increases quickly to 
a considerably higher value. This behaviour can be ex- 
plained as follows: (I) at low H concentrations diffusion 
mainly takes place in the interface region and the asso- 
ciated diffusivity is lower than in perfect bulk (red line 
in Fig. [3| due to a higher diffusion barrier in the inter- 
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FIG. 3: (Color online) Diffusion constants as a function of H 
concentrations in an idealised cubic grain for temperatures of 
(a) 300 K, (b) 600 K, and (c) 1000 K. The symbols denote 
numerical results from the kMC simulations, the dashed lines 
are obtained from Eq. ( |16[ ), respectively. The red curve rep- 
resents the diffusivity of H in perfect bcc-Fe within the dilute 
limit. 



face region; (II) when the H concentration approaches 
the number of interface sites, a dip occurs in the diffu- 
sion constant. This is due to blocking in the interface 
region, i.e. hydrogen is mainly confined to the interface 
region filling all the interface sites and effectively block- 
ing diffusion processes; (III) at high H concentrations 
(TVh > n-intf), the overall diffusivity results mainly from 
bulk diffusion, which is much faster and thus the diffusiv- 
ity increases. The dip in the diffusivity for hydrogen con- 
centrations close to the number of interface sites is less 
pronounced for higher temperatures, because at higher 
temperatures H atoms have a higher probability to es- 
cape from the interface to the bulk region. At very high 
hydrogen concentrations the diffusivity decreases again 
due to blocking in the bulk region. 

Our numerical results from the kMC simulations can 
also be described by an analytic expression. The over- 
all diffusion constant is a combination of diffusion in the 
interface and in the bulk region and can thus be approx- 
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FIG. 4: (Color online) Probability of finding hydrogen in the 
bulk region of an idealised cubic grain at T = 600 K as a 
function of H concentration. The circles represent numerical 
results from the kMC simulations, the dashed line shows the 



result from the grand canonical ensemble model, cf. Eq. (111. 



The solid blue line marks the theoretical limit ribuik/(n-buik + 

nintf). 



imated as the weighted sum over the two contributions: 



D(T,N H ) = 5> s (T)A>, s exp 



-AE S 
knT 



S S (T,A H ) 



(16) 

The sum is weighted by the probabilities, p s , to be in 
the bulk or interface region, cf. Eq. (11). In addition 



the blocking^ in the two regions is accounted for by the 
factor 



Z S (T,N H )= 1- 



iV H , s (T,iVH) 



(17) 



where N s (T,Nn) = p s (T)Nn is the average number of 
hydrogen atoms at site type s and n s is the number of 
sites of type s. Diffusion in between the two site types 
is assumed to be in equilibrium and these processes do 
not significantly contribute to the overall diffusivity. The 
prefactor -Do.s contains the geometric prefactor r s , the 
attempt frequency vq and the hopping distance ao, with 
z/ = 10 13 s _1 and a — 1.001 A. The geometric prefac- 
tors for the bulk and interface regions, Tbuik and r; nt f are 
obtained by fitting to the numerical kMC data. 

Before comparing the analytic and numerical results of 
the diffusion constants we verify that the hydrogen con- 
centration in the two regions is indeed in equilibrium. For 
this, we evaluate the probability to be in the bulk region, 
Pbuik, within our kMC simulation according to Eq. (16| 
and compare it to the analytic result in Eq. ( 11 ). In FigTR] 
the results are shown for a temperature of T = 600 K. 
The analytic result agrees well with the kMC results. 
For high H concentrations, pbuik converges to the the- 
oretical limit «.buik/("-buik + "intf)- Comparing Figs. [3] 
and [4] indicates that for a wide range of concentrations, 
hydrogen is indeed confined to the interface region, even 
at T — 600 K. Only when the interface sites are filled 
with hydrogen diffusion in the bulk region enhances the 
overall diffusivity. 



Using our kMC data at T = 300, 600, and 1000 K, the 
fitted geometric prefactors are rj nt f = 0.377 and Tbuik = 
0.644. As shown in Fig. [3] the analytic results are in 
excellent agreement with the kMC results. The fitted 
values for the geometric prefactors also agree well with 
their theoretical values. In the bulk region each site has 
four nearest neighbours yielding Tbuik = 4/6 = 0.667, 
whereas the interface region is a 3D network of 2D grain 
boundary plains with only two nearest neighbours for 
each interface site, yielding r; nt f = 2/6 = 0.333. The 
small deviation of the numerical fitted values from the 
ideal theoretical ones is mainly due to diffusion processes 
in the vicinity of neighbouring bulk and interface sites 
where the local connectivity differs from the ideal one. 



IV. IDEALISED CUBIC GRAIN WITH 
ISOLATED POINT DEFECTS 

In addition to the interface and bulk sites we next in- 
cluded point defects in the bulk region as a third site 
type into our model. Point defects such as vacancies or 
substitutional atoms are typically present in materials 
and may as well influence the diffusion of hydrogen. Our 
idealised cubic grain model with point defects contains 
400 point defects within the bulk region, 5952 interface 
sites and 42800 bulk sites. The point defect sites have 
the same geometry as bulk and interface sites, but dif- 
ferent energetics. Since hydrogen binds even stronger 
to vacancies than to interface site^U, the point defect 
sites are also considered to be more stable in our kMC 
model. The corresponding microscopic diffusion barriers 
are summarised in Table H] 

Since the point defects are not connected to each other 
there is no direct hop from one point defect to the next. 
Hence, the geometric prefactor, r trap for point defects in 
Eq. (16) is expected to be zero. However, the point de- 



fects alter the distribution of hydrogen within the system 
and thus the corresponding probabilities, p s . 

In Fig. [5] the results for H diffusivities in an idealised 
cubic grain with and without point defects are compared 
for a temperature of T = 600 K. At low concentrations 
most hydrogen atoms are trapped at the point defects 
and hence do not directly contribute to the overall dif- 
fusivity. Therefore, the diffusion constants are lower as 
compared to those in the same grain without point de- 
fects. At higher H concentrations the effect becomes less 
pronounced and eventually all point defects and interface 
sites are filled, and hydrogen atoms in the bulk region 
dominate the overall diffusion. 

For the grain model including point defects the fitted 
geometric prefactors for the interface region, bulk region 
and point defects are F; nt f = 0.376, Tbuik = 0.634, and 
0.0, respectively. The analytical results obtained 



trap 



from Eq. (16) are in excellent agreement with the kMC 
results as shown in Fig. [5] Note, that the fitted geo- 
metric prefactors for the model without point defects are 
Tintf = 0.377 and Tbuik = 0.644. The geometric prefactor 
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FIG. 5: (Color online) Diffusion constants as a function of 
H concentration at T = 600 K in an idealised cubic grain 
(blue (dark) triangles) and in and idealised grain with ad- 
ditional 400 point defects in the bulk region (orange (light), 
squares). The symbols denote results from the kMC simula- 
tions, dashed/dotted lines indicate analytic results according 
to Eq. | |16| . The red curve represents the diffusivity of H in 
perfect bcc-Fe within the dilute limit. 



FIG. 6: (Colour online) Model of an idealised layered struc- 
ture of parallel grain boundary regions in bcc-Fe. The spheres 
represent tetrahedral interstitial sites in bcc-Fe, each having 
four nearest neighbours at the same distance. Interface sites, 
intermediate sites, and bulk sites are shown in blue, red, and 
green, respectively. 



As can be seen from Eq. (21) the effect of point defects 



becomes negligible for small numbers of trapping sites, a 
small binding energy or high temperatures. 



IDEALISED LAYERED STRUCTURE 



for the bulk region decreases slightly in the presence of 
point defects, since the point defects reduce the effective 
number of equivalent nearest-neighbour sites within the 
bulk region. 

In the limit of low hydrogen concentrations and low 
temperatures Eq. ( 16 1 reduces to an expression previ- 
ously used to describe H diffusion in the presence of point 
traps in iron material J 2 * 33 l Assuming there are only bulk 
and isolated point defect trapping sites with r trap = 0.0, 
Eq. ( 16 ) reads: 



D(T, N H ) = 

Pbuik A),buik exp 



(18) 



— A£buri 

knT 



with 



Pbulk 



(19) 



"bulk 

l+cxp((_E bulk - M )/fc B T) 



l+cxp((£,„p-M)/iBr) t i + cxp((E bulk - M )/fc B T) 



For both low H concentrations and low temperatures we 
can approximate exp(_Ei — fi/k^T) 3> 1 and 



Pbulk 



"bulk 

cxp((_E b ui k -M)/fcBT) 



cxp((_E tra 



-fj.)/k B T) 
1 



cxp((_E bulk -^)/fc B T) 



■ exp 



/ A£ bind \ ' 
V fc «T ) 



(20) 



where AE 



bind 



bulk 



Strap is the binding energy of 
hydrogen to the trapping site. Furthermore, f or lo w H 
conc entr ations Sbuik — > 1 and substituting Eq. ( 20 ) into 
Eq. ( 18 ) leads to the known expressiorpl22l 



D{T) = 

A),buik exp 



k P ,T 



1 + exp 



/ AiJbind 

\ k B T J ribuik 



^trap 



(21) 
-r 



The idealised cubic grains discussed in Sections |III| 



and IV are small compared to typical grain sizes in mate- 
rials. In atomistic simulations grain boundary structures 
are thus often modelled as parallel arrangements of inter- 
face planes to describe their 2-dimensional nature. As a 
second microstructure model we investigate the diffusion 
of hydrogen within such an idealised layered structure as 
shown in Fig. [6] 

All sites in the kMC model reflect the geometry of 
tetrahedral interstitial sites in bcc-Fe. The GB region 
is represented by three layers of interface sites indicated 
by blue spheres in Fig. [6j Additionally, the three layers 
above and below the interface region are marked as in- 
termediate sites (red spheres). The structure is highly 
anisotropic, and diffusion within the interface layers and 
perpendicular to them is expected to differ significantly. 
Within the kMC simulations cells with side lengths of 
x = y = 11.328 A, and z = n x 11.328 A (with n = 1- 5) 
are used that are periodically repeated in all three di- 
mensions. 

In a first step we investigate the diffusivity of H as 
a function of the interface-interface distance. Here, the 
intermediate sites are equivalent to bulk sites, the corre- 
sponding barriers are taken from case (I) in Table [Tl The 
results of the kMC simulations at T = 600 K are shown 
in Fig. [7j In Fig [7] (a) and (b) the diagonal components 
of the diffusion tensor parallel, D\\ = 1/2(D XX + D yy ), 
and perpendicular, D± = D zz , to the interface layers 
are depicted, respectively. Diffusion perpendicular to 
the interface layers is much slower than parallel to the 
interface, but in both directions the diffusivity exhibits 
a linear dependence on the distance between the inter- 
face layers. This is due to the probability of finding 
H atoms within the bulk region, Pbuik- According to 
Eq. (16) at low H concentration (with Sbuik ~ Sj n tf ~ 1) 
the overall diffusion constant depends linearly on Pbuik 
(with pi n tf = 1 — Pbuik) where the slope depends on the 
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FIG. 7: (Colour online) H diffusivity and H distribution as a 
function of the interface-interface distances within the layered 
structure: (a) components of the diffusion tensor parallel to 
the interface, Du — 1/2(D XX + D yy ); (b) component of the 
diffusion tensor perpendicular to the interface, D± = D zz ; 
(c) probability of finding H in the bulk region, pbuik- Sym- 
bols denote results from the kMC simulations, dashed lines 
correspond to linear fits. 



difference in diffusivity in the bulk and interface region. 
Since the diffusivity in the bulk region is higher than in 
the interface region the overall diffusivity increases lin- 
early with increasing fbuik- Pbuik again increases like- 
wise linearly with increasing interface-interface distance 
as shown in Fig. Yn (c) . Since within our current model 
the condition exp(AE'bmd/fcB'7 1 ) • nintf/^buik ^> 1 holds 
(i.e. there is a notable number of interface sites with a 



significant binding energy), Eq. (20) can further be sim- 
plified 



Pbulk 



n bulk 
nintf 



exp 



— AE hind 
k B T 



(22) 



suggesting that Pbuik is proportional to the ratio of bulk 
and interface sites, nbuik/^intf , at a given temperature 
T. As nbuik increases linearly with the interface-interface 
distance, whereas ni nt f remains constant, their ratio in- 
creases linearly and so does Pbuik- 

Diffusion perpendicular to the interface layers can be 
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FIG. 8: (Colour online) Diffusion constants as a function of 
H concentration at T = 600 K in an idealised layered struc- 
ture. The diffusivities parallel (-Dm) and perpendicular (D±) 
to the interface layers are shown in (a) and (b), respectively. 
The blue triangles denote the numerical results from the kMC 
simulations, the dashed lines are the analytic results obtained 
via Eq. (16 1. The red curve represents the diffusivity of H in 



perfect bcc-Fe within the dilute limit. 



described in an even simpler model. Since the interface 
region provides traps for H atoms with an energy gain 
of 0.512 eV and H diffuses rapidly within the bulk re- 
gion, the diffusion perpendicular to the interface can be 
approximated by H hopping between two adjacent grain 
boundary planes. Within this ID-model the hoping rate 
is k ± = i/ ± e~ ABi / kBT , and 



D± = T±a ± vj_ exp 



(23) 



Here, T± = 2/2 = 1 and a± is the interlayer distance. As 
discussed, the probability of finding H in the bulk region, 
Pbuik is proportional to the interlayer distance and thus is 
the time, r, that H spends in the bulk region, yielding the 
relation pbuik ~ a± ~ t = —. Hence, the attempt fre- 
quency v_\_ is inversely proportional to a±, and as a result, 
the diffusion constant D± is proportional to aj_ , which is 
consistent with the results shown in Fig. [7] (b) . The linear 
dependence of l/v±_ on the interlayer distance is strictly 
only valid if the bulk region is much thicker than the 
interface region. With increasing interface-interface dis- 
tance and in the limit of an infinitely thin interface layer 
the effective hopping barrier AE± converges to the bar- 
rier for escaping from the interface region AE^tf^buik- 

As a second step the dependence of the diffusivity 
within the layered structure on the H concentration is 



investigated. The results of the kMC simulations are 
shown in Fig. |forT = 600 K. The simulation cell with 
z = 22.656 A contains 128 interface sites and 1408 bulk 
sites. 

Figure |8ja) illustrates the diffusivity parallel to the in- 
terface (-Dm). Similar to the results for the cubic grain 
structure the diffusivity is constant and lower than bulk 
diffusion for low H concentrations. When the number 
of H atoms approaches the number of available inter- 
face sites blocking occurs and a dip is observed for Dy. 
For even higher concentrations the diffusivity increases 
quickly due to contributions from bulk diffusion. The 
dashed line shows the results of our analytic model ob- 
tained from Eq. (16). The fitted geometric prefactors are 
FiLintf — 0.455 and T|| bu j k = 0.638, and the numeric and 
analytic results are in very good agreement. Since the in- 
terface layers can be considered as parallel, 2-dimensional 
planes, the theoretical value for the geometric prefactor 



TABLE II: Diffusion barriers, AEi, between various sites i 
within the idealised layered structure. The three cases rep- 
resent three different models to approximate different GBs in 
bcc-Fe. 



is r^intf — i = 0-5- Within the bulk region the geomet- 
ric arrangement is unchanged, i.e. T|| bu i k = | = 0.667. 
Again the fitted geometric prefactors are close to the ideal 
theoretical values. 

Figure ^b) shows the diffusivity perpendicular to the 
interface plane, D±, as a function of H concentration. For 
low concentrations diffusion perpendicular to the inter- 
faces is much slower than parallel diffusion. This is due 
to the fact that the H atoms are largely confined within 
the interface region. At large H concentrations bulk dif- 
fusion dominates and the diffusivity becomes isotropic. 
Since diffusion within the interface layers does not signif- 
icantly contribute to D±, there is also no blocking effect 
observed. Instead, the diffusivity increases smoothly as 
the number of H atoms approaches the number of in- 
terface sites and bulk diffusion becomes dominant. The 
dashed line is obtained within our analytic model. Also 
here the agreement between numerical and analytic re- 
sults is remarkable. The fitted geometric prefactors are 
r_i_,intf = 0.0 and r ± bulk = 0.470. The zero value of 
r_i_,mtf confirms that diffusion within the interface lay- 
ers does not contribute to D±. The value of Tj^buik is 
smaller than the expected value of 0.667. The deviation 
can be explained by the fact that H diffusion in the bulk 
region perpendicular to the interfaces is interrupted by 
trapping within the interface planes effectively lowering 
the overall diffusion constant. It also indicates that our 
assumption within our analytic model, that the overall 
diffusion is a weighted sum over the diffusion within the 
different regions, does not fully apply to D±. Here, the 
overall diffusion is a combination of diffusion within the 
bulk region and trapping within the interface region. 

In a third step we investigate the influence of different 
types of GBs on the diffusivity. Our DFT calculations for 
the £5 GB in bcc-Fe^ indicate that diffusion perpendic- 
ular to the GB interface proceeds via an intermediate 
site. Furthermore, the more close-packed £3[110](112) 
GB (E3 in short) in bcc-Fe exhibits rather large diffu- 
sion barriers within in the interface region as well as per- 
pendicular to itPl We consider three different cases for 



process i 


case (I) 


case (II) 


case (III) 


bulk — > bulk 


0.088 


0.088 


0.088 


bulk — > intermediate 


0.088 


0.088 


0.088 


intermediate — > bulk 


0.088 


0.288 


0.088 


intermediate — > intermediate 


0.088 


0.200 


0.088 


intermediate — > interface 


0.088 


0.088 


0.088 


interface —¥ intermediate 


0.600 


0.400 


0.500 


interface — > interface 


0.250 


0.550 


0.550 



which the barriers are summarised in Table [TTj Case (I) 
resembles our original setup where the intermediate and 
bulk sites are equivalent and the diffusion barrier between 
interface sites is larger than in the bulk region but still 
considerably lower than the diffusion barrier out of the 
interface region. In case (II) intermediate sites are intro- 
duced where the diffusion barrier between intermediate 
sites is lower but almost comparable to the escape barrier 
into the bulk region. Case (III) resembles the situation 
within the £3 GB in bcc-Fe. Both the diffusion bar- 
rier within and perpendicular to the interface region is 
rather large. As already discussed, for a temperature of 
T = 600 K and in the dilute limit for case (I) hydrogen is 
largely confined to the interface region and almost the en- 
tire diffusion takes place in the interface region, although 
the diffusion barrier is higher than in the bulk region. 
In case (II) approximately half of all diffusion processes 
are found in the intermediate region, which only provides 
about 1/5 of the available sites. This indicates that diffu- 
sion within the intermediate region is preferred, guiding 
hydrogen along the grain boundary, thereby partially al- 
leviating the trapping effect of the energetically low lying 
interface sites. For case (III) we find that the interface 
region only acts as a trap for hydrogen. Diffusion is al- 
most entirely observed in the bulk region and there is no 
enhanced diffusivity in the intermediate region, i.e. in 
this case the GB does not influence the preferred diffu- 
sion direction of hydrogen. 

Our findings for the idealised layered structures indi- 
cate that depending on the nature of the actual GB the 
interface region might cither determine the preferred dif- 
fusion direction or trap and effectively immobilise hy- 
drogen atoms at the interfaces. In the next Section we 
extend our layered model to more closely resemble the 
structure of the £5 GB in bcc-Fe. 



VI. £5 GRAIN BOUNDARY IN bcc-Fe 

The structure of the £5 grain boundary in bcc-Fe is il- 
lustrated in Fig. [9] To setup a kMC model that describes 
the diffusion of hydrogen within the £5 GB it is necessary 
to identify all stable interstitial sites for H atoms within 
this structure as well as possible diffusion processes be- 
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FIG. 9: (Colour online) Simulation cell of the E5 grain bound- 
ary in bcc-Fe. The red spheres represent Fe atoms, the GB 
interface is indicated by the blue dashed line. All intersti- 
tial sites are characterised with respect to a reference atom 
denoted by R. 



tween these sites. Here we performed extensive DFT cal- 
culations to obtain reliable values for solution energies 
and diffusion barriers. Our results are thus not depen- 
dent on any fitted parameters but all input data are ex- 
tracted from ab-initio calculations. We have investigated 
the stability of hydrogen within various interstitial sites 
in a previous studjfl^ and found that within the E5 GB 
interstitial sites close to the interface region are energet- 
ically more favourable for H atoms than the tetrahedral 
site in bulk bcc-Fe. The DFT calculations revealed 8 dis- 
tinctive interstitial sites within the £5 GB. The relative 
solution energies of H within these sites with respect to 



the most stable site (if3) are summarised in Table III 

To identify suitable transition states for diffusion pro- 
cesses between the various inter stitial sites we employed 
the nudged ela stic b and methocP^ 3 ^ as implemented in 
the VASP cod e! 36 ' 37 1 . The computational details can be 
found in Ref. [T5] We found 12 different diffusion pro- 
cesses, the corresponding barriers extracted from the 
DFT calculations are summarised in Tabic ITVl 

Based on our very detailed DFT study the kMC model 
of the £5 GB structure is constructed by mapping the 8 
identified interstitial sites onto a lattice connected by the 
12 diffusion processes listed in Table IV The (lxlxl) 
supercell of the £5 GB shown in Fig.[9jhas the dimensions 
x = 8.98 A, y = 2.84 A, and z = 71.78 A. In addition 
to the tetrahedral interstitial sites in the bulk region, the 
relative positions of the different interstitial sites close to 
the interface region with respect to a reference atom R 
(cf. Fig. [9]) are listed in Table III The reference atom R 
sits within the GB plane. Due to the mirror symmetry of 
the GB interface, for an interstitial site at (Ax, Ay, Az) 
with Az > 0.1 A, there is an equivalent site at (Ax, Ay, 
— Az). Since the £5 GB has a base-centred orthorhombic 



TABLE III: Relative positions of various symmetry- 
inequivalent interstitial sites within the £5 grain boundary 
in bcc-Fe. The positions are given with respect to the refer- 
ence atom R in Fig. [9] Solution energies are listed relative 
to the most stable interstitial site, if3. The relative solution 
energy of a tetrahedral bulk site is 0.500 eV. 



site 


Az (A) 


Ay (A) 


Az (A) 


AE (eV) 


ifl 


7.041 


0.000 


0.001 


0.050 


if2 


3.917 


0.000 


0.007 


0.269 


if.3 


5.522 


0.000 


0.075 


0.000 


iml 


4.249 


0.000 


0.659 


0.237 


im2 


1.287 


0.000 


0.759 


0.208 


im3 


7.591 


0.000 


1.204 


0.315 


im4 


6.854 


-0.001 


1.504 


0.373 



TABLE IV: Hydrogen diffusion barriers for possible transi- 
tions between interstitial sites in the E5 grain boundary in 
bcc-Fe. Values are given for both the forward and backward 
process. The process im3 4-¥ im2 takes place between second 
nearest neighbours, while all other transitions are between 
nearest neighbour sites. 



process 


A_Ef orwa rd (eV) 


A-Ebackward (©V) 


if3 <-» ifl 


0.118 


0.068 


ifl im3 


0.339 


0.074 


im3 <-» im2 


0.109 


0.216 


im2 <->• bulk 


0.383 


0.091 


if3 im4 


0.429 


0.056 


if 2 iml 


0.020 


0.052 


iml «-> if3 


0.102 


0.339 


im4 <-¥ bulk 


0.201 


0.074 


bulk ■<-> bulk 


0.088 


0.088 


if3 if3 


0.250 


0.250 


if3 <r> im2 


0.241 


0.033 


im4 o im3 


0.030 


0.088 



structure, for an interstitial site at (Ax, Ay, Az) there 
exists an equivalent site at (Ax + x/2, Ay + y/2, Az). 
For the kMC simulations we used a model corresponding 
to a (2 x 6 x 1) supercell that is repeated periodically in 
all three dimensions. The model contains 48 ifl, 48 if2, 
48 if3, 96 iml, 96 im2, 96 im3, 96 im4, and 10224 bulk 
sites. 

The interstitial sites and microscopic diffusion pro- 
cesses create a complex diffusion network for hydrogen 
that is highly anisotropic. To illustrate this network the 
contribution of different processes to the overall diffusiv- 
ity is shown in Fig. [10] for low H concentrations. The 
relative line thickness corresponds to the relative prob- 
abilities of observing H atoms diffuse along this connec- 
tion. For all three temperatures T = 300, 600 and 1000 K 
diffusion within the bulk region is negligible. Thus dif- 
fusion in z-direction is slow since for this hydrogen has 
to leave the interface region and cross the bulk region 
towards the next interface. The plots on the left side of 
Fig. 10 show a view along the y-axis, [001]-direction, il- 



lustrating the diffusion network in x-direction. Although 
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FIG. 10: (Colour online) Diffusion network within the £5 
GB. The red spheres represent Fe atoms, and the blue lines 
illustrate the diffusion probabilities for various H diffusion 
paths. 



there is a thick connection between the if3 and ifl site 
there exists no continuous diffusion pathway within the 
interface. To diffuse in x-direction hydrogen has to leave 
the interface taking less favourable paths. Thus diffusiv- 
ity in x-direction is low for small hydrogen concentrations 
as also shown in Fig. |11| For the diffusion in y-direction 
hydrogen can continuously move along if3 sites without 
leaving the interface region as shown in the plots on the 
right side of Fig. 10 (view along x-axis, [130]-direction). 
With increasing temperature the diffusion network be- 
comes more isotropic within the interface plane. This is 
also reflected in the different components of the diffusion 

For T = 300 K at low H con- 



tensor shown in Fig. 11 
centrations D xx and D yy differ by almost two orders of 
magnitude, whereas at T = 1000 K they are comparable. 

As discussed in the previous sections, for the idealised 
interface structures the diffusivity shows a strong depen- 
dence on the hydrogen concentration. In Fig. 11 the di- 
agonal components of the diffusion constant tensor as 
a function of the hydrogen concentration are shown for 
T = 300, 600, and 1000 K for the £5 GB. For low H 



concentrations diffusion in the interface region dominates 
and the diffusivity is anisotropic as expected from the 
discussion of the diffusion network. For large H concen- 
trations diffusion is again isotropic since bulk diffusion 
dominates. The overall diffusion constants measured in 
the kMC simulations are the result of the statistical in- 
terplay of all possible diffusion processes. In our simple 
analytic model diffusion within a certain site type is as- 
sociated with one specific barrier, whereas in the £5 GB 
there are several processes and thus barriers associated 
with one site type and diffusion within the interface layer 
proceeds via several site types. It is thus not possible to 
accurately describe the diffusivity within our simple an- 
alytic model in Eq. (16). 



Nevertheless, in order to obtain an understanding of 
how much the diffusion within the GB interface and 
the bulk region contribute to the diffusion constant, we 
approximate the diffusivity within £5 GB by selected 
and/or combined processes within the analytic model. 

H diffusion along the x-direction, and within the GB 
interface follows the path if3— HmA— >im3— >im2— »if3 (cf. 
Fig. 10 ). We assume that this diffusion path is dominated 



by the step if3— >im4 with a barrier of 0.429 eV. Following 
Eq. (16 1, we approximate the diffusivity in x-direction, 



(24) 



;<ii:;iniV, .-.■„.; <-l> ( A ^' f ^ > ' m4 ) -ir.-il-T. X u i 



+Pbulk(T)D ,xx : bu\k exp ( ^y"' k J -bulk(T, A H ) , 

with Dq^xx.s — T XXtS a%. x s vq. The other contributions to 
the sum are neglected. Still, as in the case of isolated 
point defects, the other interstitial sites influence the hy- 
drogen distribution and thus the corresponding probabil- 
ities in Eq. (24). The hop lengths between two adjacent 
if3 sites along the x-direction is a X x,if3 — 4.486 A. 

Along the y-direction, the overall diffusion constant, 
Dyy, is similarly approximated within our analytical 
model. Diffusion in y-direction either takes place within 
the bulk region or in between if3 sites. For the hop be- 
tween two adjacent if3 sites along the y-direction, the hop 
length is a w ,if3 = 2.837 A and the barrier A£ , jf 3 _j.if 3 = 
0.250 eV, the other diffusion processes are not taken into 
account. 

As discussed in Sec. [Vj for H diffusion along the in- 
direction (perpendicular to the interface), the GB inter- 
face acts like a point trap. All geometric prefactors for 
diffusion between interface sites are essentially zero, and 
only bulk diffusion contributes to D zz . 

The analytic results are shown together with the kMC 
results in Fig. [IT] For D yy the agreement is excellent, 
indicating that diffusion in bulk and along if3 sites are 
indeed the dominating processes. For D xx and D zz the 
analytic results deviate from the numeric data, especially 
for lower temperatures. Still, considering the simplicity 
of our analytic model the agreement is quite remarkable 
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FIG. 11: (Colour online) Components of the diffusion con- 
stant tensor as a function of H concentration within the £5 
GB in bcc-Fe. Results are shown for (a) 300 K, (b) 600 K, 
and (c) 1000 K. The kMC results for D xx , D vy , D zz , and the 
overall diffusivity D are represented by blue squares, green di- 
amonds, orange circles, and black triangles respectively. The 
analytic results are shown as dashed lines in the correspond- 
ing colour. The red curve represents the diffusivity of H in 
perfect bcc-Fe within the dilute limit. 



region is actually smaller than in between interface sites, 
i.e. -D TO ,buik < Dyy^fs with £>„„ = v,fzD, 
cf. 



yy — Pif3-^yy,if3 ^if3 

p4|. The diffusion barrier 



Pbulk -Dyy, bulk Sbulk, cf. Eq 

between if3 sites is larger than between bulk sites, but at 
high temperatures the prefactor, which is about a factor 
of 10 larger for diffusion between if3 sites, dominates. 
Thus at high concentrations where bulk diffusion has 
the largest contribution the overall diffusivity, D yy , de- 
creases. At all temperatures, diffusion perpendicular to 
the interface planes, D zz , is slowest, due to trapping of H 
within the grain boundary region. With respect to diffu- 
sion in x-direction (D xx ) the kMC results for T = 300 K 
suggest that the if3 site likewise acts as point trap. At 
this low temperature H atoms that occupy if3 sites are 
nearly immobile in x-direction. Once the H concentra- 
tion reaches that of if3 sites, diffusion in x-direction may 
occur following the path im2 ■<-> bulk — > im2 or im4 — >• 
bulk — > im4 with smaller diffusion barriers of 0.38 eV 
and 0.21 eV, respectively. These effectively lower bar- 
riers are also consistent with the increase and relatively 
large value of D xx with increasing H concentration as 
shown in Fig. 
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(a). For high temperatures and high 
H concentrations bulk diffusion dominates and the diag- 
onal components of the diffusion tensor are equivalent, 
indicating an isotropic diffusion behaviour. In all cases 
the overall diffusion in the grain boundary structure is 
slower than diffusion in perfect bcc-Fe bulk. 

It is apparent that for the more realistic model of the 
E5 GB in bcc-Fe our analytic model is too simple to 
fully describe the diffusion of hydrogen. Nevertheless, 
the observed trends at high temperatures may be under- 
stood qualitatively with a simple analytic approximation. 
Analysing the kMC results a diffusion network was es- 
tablished and the processes that dominate the diffusion 
under various conditions were extracted yielding a more 
detailed understanding of hydrogen diffusion within grain 
boundary structures. 



for T = 600 and 1000 K. 

The fitted values of the geometric prefactors are in in- 
direction r^ifs = 1.531, r xX; bulk = 0.747, in y-direction 
ryy,if3 = 1-074, Tyy^uik = 0.701, and in z-direction 
r z z,buik = 0.705. All bulk values are relatively close to 
the theoretical value of 0.667. Diffusion in y-direction 
along if3 sites corresponds essentially to a ID chain with 
two nearest neighbours, yielding a theoretical value of 
^yy,ia — 2/2 = 1. The good agreement between fitted 
and theoretical values for the geometric prefactors in y- 
direction also indicates that the simple analytic model is 
suitable here. 

Diffusion within y-direction, D yy , exhibits a small 
blocking effect as seen within our idealised models, which 
is again most pronounced for low temperatures and van- 
ishes at higher temperatures. At T — 1000 K there 
is no dip, but D yy somewhat decreases with increas- 
ing H concentration. This is due to the fact that at 
1000 K the effective diffusion constant within the bulk 



VII. CONCLUSIONS 

Employing kinetic Monte Carlo simulations we have 
studied hydrogen diffusion within various models repre- 
senting different arrangements of grain boundaries and 
point defects in bcc-Fe. The defect regions exhibit inter- 
stitial sites with a significantly lower solution energy for 
H atoms, effectively acting as trapping sites. Within an 
idealised cubic grain structure we observe a characteristic 
behaviour of the diffusion tensor as a function of hydro- 
gen concentration. At low concentrations H is confined 
to the interface region and the diffusivity is low as com- 
pared to diffusion in perfect bcc-Fe bulk. As the number 
of H atoms approaches the number of interface sites the 
diffusion constant drops due to blocking of available in- 
terstitial sites. At large H concentrations bulk diffusion 
dominates the behaviour and a significant increase of the 
diffusivity is observed. Additional point defects lower the 
diffusivity for small H concentrations, but do not change 
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the behaviour for larger concentrations. 

Within a layered arrangement of grain boundary 
planes the diffusion is anisotropic. Parallel to the in- 
terface diffusion is similar to the one observed within the 
grain structure. Perpendicular to the interface diffusion 
is much slower and can effectively be described by a ID 
model of H atoms hopping between neighbouring inter- 
face planes. The effect of the grain boundary on the 
diffusion of hydrogen strongly depends on the actual so- 
lution energies and diffusion barriers between different 
interstitial sites. 

The more detailed model of the £5 GB in bcc-Fe 
showed that the overall diffusion is a complex interplay 
of various microscopic diffusion processes. Depending 
on the conditions (hydrogen concentration, temperature, 
diffusion direction) different processes dominate the dif- 
fusion resulting in a complex diffusion network. Still, the 
general trends with respect to temperature and hydro- 
gen concentration may be understood from a simplified 
analytic model. 

We have derived a simple, analytical expression for hy- 
drogen diffusion within microstructures that consist of 
several distinctive regions (such as bulk, interfaces, point 
defects). The analytic model is in very good agreement 



with the numerical results for the idealised structures. 
For the more detailed model of the S5 GB the analytic 
model only works for conditions where the diffusion is 
dominated by a few, specific processes, but naturally it 
fails to capture the more complex interplay between a 
number of different microscopic diffusion processes. 

In all structures the diffusivity is lower than in perfect 
bcc-Fe bulk, indicating that the grain boundary regions 
do not serve as fast diffusion channels. Nevertheless, at 
low concentrations hydrogen is confined to the interface 
region, i.e. the arrangement of grain boundary planes 
and thus the microstructure significantly influences the 
preferred diffusion direction. 
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